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Abstract: We present a nonlocal electrostatic formulation of nonuniform ions and water molecules with in¬ 
terstitial voids that uses a Fermi-like distribution to account for steric and correlation effects in electrolyte 
solutions. The formulation is based on the volume exclusion of hard spheres leading to a steric potential and 
Maxwell’s displacement field with Yukawa-type interactions resulting in a nonlocal electric potential. The 
classical Poisson-Boltzmann model fails to describe steric and correlation effects important in a variety of 
chemical and biological systems, especially in high field or large concentration conditions found in and near 
binding sites, ion channels, and electrodes. Steric effects and correlations are apparent when we compare 
nonlocal Poisson-Fermi results to Poisson-Boltzmann calculations in electric double layer and to experimen¬ 
tal measurements on the selectivity of potassium channels for I< + over Na + . 
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1 Poisson-Fermi Theory 

Continuum electrostatic theory is a fundamental tool for studying physical and chemical properties of elec¬ 
trolyte solutions in a wide range of applications in electrochemistry, biophysics, colloid science, and nanoflu¬ 
idics [1-8]. For over a century, a great deal of effort has been devoted to improving the Poisson-Boltzmann 
(PB) theory of continuum models for a proper description of steric (or finite size) and correlation (or nonlocal 
screening, polarization) effects in electrolytes [9-24]. We present a continuum model with Fermi-like distri¬ 
butions and global electrostatic screening of nonuniform ions and water molecules to describe the steric and 
correlation effects, respectively, in aqueous electrolyte solutions. 

For an electrolytic system with K species of ions, the entropy model proposed in [20] treats all ions and 
water of any diameter as nonuniform hard spheres and regards the water as the (K + l) th species. It then 
includes the voids between these hard spheres as the (K + 2) th species so that the total volume V of the 
system can be calculated exactly by the identity 

K+l 

V = ^2 V i N i + V K+2> (1) 

2=1 

where V K+2 denotes the total volume of all the voids, v,- = 47ra?/3 that gives the volume of each sphere 
with radius a,-, and JV, is the total number of the i th species particles. In the bulk solution, we have the bulk 
concentrations Cf = ^ and the bulk volume fraction of voids r B = . Dividing the volume identity (1) 
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by V, r B = 1 - is expressed in terms of nonuniform v, and Cf for all particle species. If the sys¬ 

tem is spatially inhomogeneous with variable electric or steric fields, as in most biological and technological 
systems, the bulk concentrations then change to concentration functions C,(r) that vary with positions, and 
differ from their constant values Cf at location r in the solvent domain Q s . Consequently, the void volume 
fraction becomes a function T(r) = 1 - )T^ 1 v,Ci(r) as well. 

For an electrolyte in contact with electrodes or containing a charged protein, an electric field E(r) in the 
solvent domain Q s is generated by the electrodes, ionic (free) charges with a displacement field D(r), and 
bound charges of polar water with a polarization field P(r). In Maxwell’s theory, these fields form a constitu¬ 
tive relation 

D(r) = e 0 E(r) + P(r) (2) 

that yields the Maxwell’s equation V • D(r) = pj(r) = \ <?/C;(r), Vr e D s , where e 0 is the vacuum permit¬ 

tivity and q t is the charge on each i species ion [25, 26]. The electric field E(r) is thus screened by water (in 
what might be called Bjerrum screening) and ions (in Debye screening) in a correlated manner that is usually 
characterized by a correlation (screening) length A [4,12,27]. The screened force between two charges in ionic 
solutions (at r and r in D s ) has been studied extensively in classical field theory and is often described by the 

, -|r-r'|M 

screening kernel G(r - r) = [1], which is called a Yukawa-type kernel in [4, 24, 28, 29], and satisfies 

the partial differential equation (PDE) [24,28,29] 

-2\G(r - r ) + G(r - r ) = 5(r - r ), r, r e R 3 (3) 

A z 

in the whole space R 3 , where A = V • V is the Laplace operator with respect to r and <5(r - r ) is the Dirac delta 
function at r . The potential <p( r) defined in D(r) = -e s £oV0(r) thus describes a local potential of free ions 
according to the Poisson equation 

-e s e 0 V • V0(r) = p,(r), Vr e Q s , (4) 

where e s is a dielectric constant in the solvent domain. This local potential does not deal with the long range 
correlations between different ions or between ions and polar water in high field or crowded conditions under 
which the size and valence of ions and the polarization of water play significant roles [4,7,12,18-21,27]. We 
introduce a global electric potential 0(r) of the screened electric field E(r) to deal with the correlation and 
polarization effects in electrolyte solutions. It is written as a convolution of the local potential (p(r') with the 
global screening kernel G(r - r ), i.e., 

0(r) = J (5) 

However, it would be too expensive to calculate (p(r) using this equation. Multiplying Eq. (3) by <p(r) and then 
integrating over R 3 with respect to r [24, 28, 29], we obtain 

-A 2 Acp{ r) + 0(r) = 0(r), r e O s , (6) 

a PDE that approximates Eq. (5) in a sufficiently large domain Q s with boundary conditions <p( r) = <p(r) = 0 
on dO s and describes the relation between global (p(r) and local <p{ r) electric potentials. From Eqs. (4) and 
(6), we obtain the fourth-order PDE 

e s e 0 A 2 A{A(p(r)) - e s e 0 A(p(r) = p^r), r e D s . (7) 

Thus, when we set E(r) = -V$(r), we can use Eq. (2) to find the polarization field 

P(r) = £s£ 0 A 2 V{Acp{ r)) - (e s - l)e o V0(r). (8) 

If A = 0, we recover the standard Poisson equation (4) and the approximate polarization P = e 0 (e s - 1)E with 
the electric susceptibility e s -1 (and thus the dielectric constant e s ) if water is treated as a time independent. 
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isotropic, and linear dielectric medium [26]. In this case, the field relation D = e s £oE with the scalar constant 
permittivity e s £o is an approximation of the exact relation Eq. (2) due to the simplification of the dielectric 
responses of the medium material to the electric field E. 

We introduce a Gibbs free energy functional for the system as 


F( C, (p ) = F eI ( C, cp ) + Fen(C), (9) 

Feii C, (p) = J f P\(pdr =jj p^pxdr, 

Os O s 


Fen( C) 






dr, 


where F e ;(C, (p) is an electrostatic functional, F en (C) is an entropy functional, C = (Ci(r), C2W, •••, C/c+i(r)), 
Vo = ( vJiY v,) /(K + 1) an average volume, and L~ l is the inverse of the self-adjoint positive linear operator 
L = e s e 0 X 2 AA - e s e 0 d [24, 28,29]. Taking the variations of F(C, (p ) at cp gives Eq. (7). Taking the variations of 
F(C, <p ) at C;(r) [24,28,29] yields Fermi-like distributions 


C ; (r) = Cf exp (~Pt(p(r) + ^S trc (r) j , 


S trc (r) = In 



( 10 ) 


for all i = 1, • • • , K +1 (ions and water), where /?,■ = qJk B T with q K+1 = 0, k B is the Boltzmann constant, and 
T is an absolute temperature. The distribution is of Fermi type since it saturates. All concentration functions 
C,(r) < i [20], i.e., C,(r) cannot exceed the maximum value 1/v,- for any arbitrary (or even infinite) potential 
(p(r) in the domain Q s - In these Fermi distributions, it is impossible for a particle volume v,- to be completely 
filled with particles, i.e., it is impossible to have v,C,(r) = 1 (and thus F(r) = 0) since that would yield S trc (r) = 
-00 and hence v ; - C ; (r) = 0, a contradiction. For this reason, we must include the void as a separate species 
if water and ions are all treated as hard spheres [20], Here we do represent water and ions as spheres. Our 
approach allows other shapes to be used as well but that is not done here. 

The nonlocal Poisson-Fermi (PF) Eqs. (7) and (10) have new features, of some importance. 

(i) The Fermi-like distribution of uniform spherical ions with voids in ionic liquids was first derived in 
[6,12] using Bikerman’s lattice model [30]. The entropy functional in [12] involves a reciprocal of a minimum 
volume v with a volume fraction & that is an empirical fitting parameter and may have to be unrealistically 
large to fit experimental data in some applications [6]. It is shown in [20] that the entropy functional in [12] 
does not directly yield classical Boltzmann distributions C,(r) = Cf exp (—/?,0 (r) ) as v — > 0. It can be easily 
seen from (10) that the entropy functional F en { C) in Eq. (8) yields Boltzmann distributions as V; — > 0 for all 
i. Our derivation of F en (C) does not employ any lattice models but simply uses the volume equation (1). The 
functional F en (C) is a new modification of that in [20], where the classical Gibbs entropy is now generalized 
to include all species (ions, water, and voids) in electrolytes with the same entropy form. In fact, our F(r) is an 
analytical extension of the void fraction 1 - O in Bikerman’s excess chemical potential [6], where all volume 
parameters v,- (including the bulk fraction F B ) are physical not empirical. The adjustable parameter in our 
model is the correlation length A = 2a,- depending on the ionic species i of interest [17,20]. The PF model was 
first proposed in [6] without derivation and has been shown to produce results that are not only comparable 
to molecular dynamics (MD) simulation or experimental data but also provide insight into nonlinear proper¬ 
ties of concentrated electrolytes and ionic liquids [31-34]. Here, we formally derive the PF model for general 
electrolytes using a hard-sphere instead of lattice model with the steric potential S trc (r) first introduced in [17]. 
As compared with lattice models in [6] and demonstrated, for example, in [18, 20, 35], hard-sphere models 
significantly improve the agreement between simulation and experiment. 

(ii) The fourth-order PDE (7) is similar to those in [12, 27] used in previous papers [17, 20]. However, the 
physical origin of the PDE is different. In [27], the global convolution is performed only on the charge density 
of point-like counter ions in contrast the convolution Eq. (5) to the potential (p(r) by Eq. (4) that is generated 
by all spherical ions. In [12], a derivation for electrolytes or ionic liquids with steric effects is given from a 
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free energy function of a gradient expansion of nonlocal electrostatic energy in terms of A(p. Here, the fourth- 
order PDE is derived directly from Maxwell’s equation with the Yukawa screening kernel. Our result does not 
depend on the convergence properties of an expansion of nonlocal electrostatic energy. The use of Green’s 
function as a kernel in all previous works and the present work is the same although the convoluted function 
is quite different. From the definition of the displacement, electric, and polarization fields in Eq. (2), the 
association of the displacement field and free charges as in Eq. (4) [25, 26], and our derivation leading to the 
polarization field in Eq. (8), the new convolution formula (5) is more consistent with Maxwell’s theory since 
our derivation is more concise without using the decomposition of long and short Coulomb interactions [27] or 
the gradient expansion of non-local electrostatic energy [12]. The decomposition approach requires a mean- 
field approximation and a strong-coupling expansion for the long-distance and short-distance interactions, 
respectively [27]. The fourth-order PDE is only the first order approximation of the energy expansion [12], 

(iii) Eq. (7) defines a dielectric operator e = e s e 0 (l - A 2 A) that in turn implicitly yields a dielectric func¬ 
tion e(r) as an output of solving Eq. (7) [12, 20]. A physical interpretation of the operator was first introduced 
in [12] to describe the nonlocal permittivity in a correlated ionic liquid. The exact value of e(r) at any r £ 
D s cannot be obtained from Eq. (7) but can be approximated by the simple formula e(r) = e 0 + C]| 2 o(r)(e s - 
l)e 0 /Cg i0 since the water density function Cn 2 o(r) = C/<-+ i(r) is an output of Eq. (10). This formula is only for 
visualizing (approximately) the profile of e or e. It is not an input of calculation. The input is the operator 
? = e s £o (l - A 1 A) (with its dependence on the input parameter the correlation length A). Therefore, the PF 
Eq. (7) accounts for electrostatic (ionic charges in Eq. (4)), correlation (A in Eq. (5) and in e), polarization (Eq. 
(8)), nonlocal (Eq. (5)), and excluded volume (Eq. (1)) effects in electrolytes with Yukawa shielding with only 
one parameter A. 

(iv) The factor v,/ v o multiplying S trc (r) in Eq. (10) is a modification of the unity used in our previous 

work [20]. The steric energy -TiS trc (r)/c B T [20] of a type i particle depends not only on the emptiness (T(r) = 
1 - v,-C,-(r)) (or equivalently crowding) at r but also on the volume v,- of each type of particle. If all v,- are 

equal (and thus v, = v 0 ), then all particle species at any location r £ D s have the same steric energy and the 
uniform particles are indistinguishable in steric energy. The steric potential is a mean-field approximation of 
Lennard-Jones (L-J) potentials that describe local variations of L-J distances (and thus empty voids) between 
every pair of particles. L-J potentials are highly oscillatory and extremely expensive and unstable to compute 
numerically. This dependence is worse than a nuisance. The L-J potentials are not well defined experimentally 
(e.g., the combining rules are not well defined, whether Lorentz-Berthelot or Kong rules are used). They are 
likely to depend on ionic species, concentration and perhaps other variables. It is dangerous to have a model 
that depends sensitively on parameters that are unknown experimentally. 

(v) The global convolution in Eq. (5) may seem similar to those in [4, 24, 28, 29] but it is not. The Poisson 
equation (4) takes the place of the Fourier-Lorentzian (FL) equation — an integro-differential equation — in 
the previous work [4, 24, 28, 29] in which the local potential <p(r ) in Eq. (5) is replaced by a global electric 
potential <J>(r). Moreover, the factor 1/A 2 in Eq. (5) is replaced by (e s - eJ/A 2 , where £« and A are both 
adjustable parameters. The choice of three parameters e s , £<*,, and A in the FL model is reduced to only one A 
here. 


2 Results 

The nonlocal PF model is first compared with the modified PB model (mPB) of Borukhov et al. [36] in which 
ions are treated as cubes without considering void and correlation effects. The classical PB model (with A = 
v,- = 0 for all i, i.e., no size, void, and correlation effects) produces unphysically high concentrations of anions 
near the charged wall atx = 0 as shown by the dashed curve in Fig. 1. The surface charge density is le/(50A 2 ) 
in contact with a 0.1 M C 4 A aqueous electrolyte, where the radius of both cations and anions is a = 4.65 Ain 
contrast to an edge length of 7.5 A of cubical ions in [36], e is the proton charge, and e s = 80. The multivalent 
ions A 4- represent large polyanions adsorbed onto a charged Langmuir monolayer in experiments [36]. The 
dotted curve in Fig. 1 is similar to that of mPB in [36] and was obtained by the PF model with the size effect 
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but without voids and correlations, i.e., A = 0, V K+2 = 0 (no voids), and v K+1 = v H2 o = 0 (water is volumeless 
and hence r B = 1 - J^f= i wCf' s the bulk water volume fraction). The voids ( V K+2 f 0) and water molecules 
(V[| 2 o f 0) have slight effects on A 4- concentration (because of saturation) and electric potential (because 
water and voids have no charges) profiles as shown by the thin solid curves in Figs. 1 and 2, respectively, when 
compared with the dotted curves. However, correlations (with A = 1.6a [12]) of ions have significant effects 
on ion distributions as shown by the thick solid and dash-dotted curves in Fig. 1, where the Stern like layer 
is on the order of ionic radius a [37] and the overscreening layer [12] (C A /, (x) » 0) of excess coions (Cc> (x) > 
C B + = 0.4 M) is about 18 A in thickness. The Stern like boundary layer is an output (not a prescribed condition) 
of our model and has no special significance here, except historical. Boundary layers are found frequently 
where physical properties change discontinuously or boundary conditions are imposed by mathematicians 
or physical systems. The electric potentials 0(0) = 5.6 at x = 0 and 0(11.5) = -1.97 k B T/e in Fig. 2 obtained 
by PF with voids and correlations deviate dramatically from those by previous models for nearly 100% at 
x = 0 (in the Stern layer) and 70% at x = 11.5 A (in the screening layer) when compared with the maximum 
potential 0(0) = 2.82 k B T/e of previous models. The PF potential depth 0(11.5) = -1.97 k B T/e of the over¬ 
screening layer is very sensitive to the value of the correlation parameter A = 1.6a since the depth tends to 
zero as a —> 0. The parameter value can only be justified with experimental or molecular dynamics (MD) data 
[17,18, 20,21]. 

The computational domain j(x, y, z) : 0 < x < 40, - 5 < y < 5, - 5 < z < 5 A j for the results shown in 
Figs. 1 and 2 was chosen heuristically so that it is large enough to observe a zero Dirichlet boundary condition 
for 0(r) at x = 40 A within the accuracy to the third decimal place for five grid points near and at x = 40 A for 
all potential profiles in Fig. 2. Zero Neumann boundary condition is given on the four walls adjacent to the 
charged wall (x = 0) on which the condition is -e s V0 • n = a with n = (-1,0,0) and a = le/(50A 2 ). The 
resulting water density (not shown) from Eq. (10) is approximately equal to the bulk density C \ 20 = 55.5 M 
at those grid points. 

In implementation, the fourth-order PF equation is reduced to two second-order PDEs in [17] where 
boundary conditions for those equations can also be found. The computational efficiency of the PF model 
for solving linear systems is thus commensurable to that of standard or modified PB models. However, since 
the steric potential S trc (r) in Eq. (10) depends sensitively on the unknown electric potential 0(r) in a highly 
nonlinear manner and is obtained by solving the PF equation (7) coupled with Eq. (10) using exact or inexact 
Newton’s iterative method [17], the computational cost of the PF model for nonlinear solver to achieve self- 
consistent convergence is in general much higher than that of PB models. For example, it took 42s of the CPU 
time for the PB profile (dashed curve) in Fig. 2 comparing with lh 22min 50s (about 118 times costly) for the 
PF profile (thick solid curve) on a laptop with Intel Core i5-3230M CPU. The mesh size of the FD grid for this 
example is 0.25 A that yields a size of 270,641 for the matrix system. The error tolerance was set to 0.0001 and 
0.005 for the linear and nonlinear solvers, respectively. It was very hard to converge for the PF model because 
the surface charge density a is very high resulting in the saturation (i.e., the void fraction function Hr) in Eq. 
(10) is very close zero) of ionic concentrations as shown in Fig. 1 and the ion is very large resulting in strong 
nonlinearity in terms of the correlation parameter A. Obviously, there is much room to improve the numerical 
method of the nonlinear solver used here and proposed in [17]. 

We next show the size effect of different ions with voids in a biological system. One of the main uses of 
ionic solution theory is the understanding of ions in protein channels that control a wide range of biological 
functions [7, 38, 39]. We consider the crystal structure of the potassium channel KcsA (PDB [40] ID 3F5W 
[41]) as shown in Fig. 3, where the spheres denote five specific cation binding sites (SO to S4) [42]. The crystal 
structure with a total of JV = 31268 charged atoms is embedded in the protein domain O p while the binding 
sites are in the solvent domain D s . The exquisite selectivity of K + (with a K + = 1.33 A) over other cations such 
as Na + (a Na + = 0.95 A) by potassium channels is an intriguing quest in channology. It can be quantified 
by the free energy (G) difference of I< + and Na + in the pore and in the bulk solution [42-45]. Experimental 
measurements [43-45] showed that the relative free energy [42] 

AAG( K + -> Na + ) = [G p0 re(Na + ) - G bulk (Na + )] - [G p0 re(K + ) - G bulk (K + )] (11) 
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Figure 1: Concentration profiles of anions C A 4- (x) and cations Cc+(x) obtained by various models in a C 4 A electrolyte solution 
with a positively charged surface at x = 0. 



Figure 2: Electric potential profiles <p{x). 


is greater than zero in the order of 5-6 kcal/mol unfavorable for Na + . The electric and steric potentials at the 
binding site S2 (as considered in [42]) can be calculated on the atomic scale using the following algebraic 
formulas 


$S2 


4 ne 0 


, b fV 

EE 


<7; 


6 e p (r)\c } -A k \ 


<?S2 

£h«S2 


otrc 
> *- > S2 


In 


1 


_ VS2 
^S2 

fB 9 


( 12 ) 


where S2 = Na + or K + (the site is occupied by a Na + or a K + ), qj is the charge on the atom; in the protein given 
by PDB2PQR [46], e p {r ) = 1 + 77r/(27.7 + r) [47], r = |c ; — Cs 2 |» Cj is the center of atom;, A k is one of six 
symmetric surface points on the spherical S2, e b = 3.6, and P S2 = 1.5v K + is a volume containing the ion at 
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S2. The crucial parameter in (12) is the ionic radius a S2 = 0.95 or 1.33 A (also in |Cj - A k \) that affects (p s2 
very strongly but SgJ weakly. We obtained AAG = 5.26 kcal/mol in accord with the MD result 5.3 kcal/mol 
[42], where G p0 re(Na + ) = 4.4, G bulk (Na + ) = -0.26, G p0 re(K + ) = -0.87, G bulk (I< + ) = -0.27 kcal/mol, </> Na+ =7.5 



Figure 3: The crystal structure of the K + channel KcsA (PDB ID 3F5W) [41] with five cation binding sites SO, SI, S2, S3, and S4 
[42] marked by spheres. 


For the above results, the bulk energies G bu i k (Na + ) and G bu i k (K + ) were obtained from the experimental 
results in [48]. The pore energies G por e(Na + ) and G pore (K + ) were obtained by the formula G pore (S2) = qs 2 ( Ps 2 ~ 
^S^k B T, where S2 denotes Na + or K + and T = 298.15. The values of cp S2 (he. </> N a + and tp K <) and Sgf (SjiJ a , 
and S^T) in G p0 re(S2) were obtained using the equations in (12) with a$ 2 = 0.95 (Na + ) or 1.33 A (K + ). The 
electric potential (p S2 is generated by the charges qj of all N = 31268 atoms in the protein and the charge q S2 
in the bound ion Na + or K + at S2. The steric potential Si/) depends on the volume of the bound ion Na + or K + 
and the assumed volume V$ 2 of the binding site, which is adjustable. However, a slight change of Vs 2 does 
not affect significantly the steric energy difference between Na + and K + as shown by the steric formula in (12). 
Therefore, the change does not alter significantly the selectivity result obtained from Eq. (11). Other values 
(not shown) of cp si at other sites i = 0,1,3, and 4 can be calculated in the same way. These </> Si values in the filter 
domain Qf containing the five binding sites in Fig. 3 can be used as Dirichlet boundary data for the PF model 
(7) in the rest of the solvent domain 12 s /£3y [49] to obtain a continuous electric potential function <p(i) (not 
shown) in the whole domain D = O s u D p , where Q p denotes the domain of the protein and the membrane. 
Note that another Poisson equation is defined in Q p and some conditions for <p( r) should be imposed on 
the interface dQ p between Q s and Q p . We refer to [17, 49] for more details on the model formulation and 
numerical methods. 
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3 Conclusion 

In summary, a nonlocal Poisson-Fermi model is proposed to describe global electrostatic and steric effects 
that play a significant role of ionic activities in electrolyte solutions especially in high field or large concen¬ 
tration conditions. The model is based on Maxwell’s field theory and nonuniform hard spheres of all ions 
and water molecules with interstitial voids. The Fermi-like distribution formula can describe the distribution 
of nonuniform spherical ions and water molecules with interstitial voids. The steric potential is a mean-field 
description of Lennard-Jones potentials between particles. Poisson’s equation is self-consistent with Fermi 
distributions and global electrostatics. The present theory can be used to describe complex functions of bio¬ 
logical or chemical structures on both atomic and macroscopic scales. Comparisons with experimental data 
are promising but incomplete. 
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